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1.  Introduction 


Our  work  has  two  main  thrusts. 

The  first  thrust  is  to  develop  high  quality  numerical  tools  for  resolu¬ 
tion  of  complex  (often  stochastic  and  often  dynamic)  geometries  within 
a  simulation.  High  order  geometry  has  taken  our  front  tracking  code 
FronTier  to  a  new  scientihc  level.  Novel  theories  of  stochastic  conver¬ 
gence  have  been  applied  to  turbulent  combustion.  A  new  and  very 
promising  method  for  modeling  brittle  fracture  was  developed 

The  second  thrust  is  to  explore  the  use  of  these  tools  in  a  series  of  vital, 
ARO  relevant  applications:  brittle  fracture,  parachute  drop,  turbulent 
combustion  and  robust  hnite  element  methods  for  unstructured  grids. 

2.  Numerical  Tools  for  Complex,  Dynamic  Geometry 

2.1.  Weighted  Least  Squares  over  Discrete  Surfaces.  During 
this  project,  we  have  developed  a  unified  theoretical  framework  based 
on  the  weighted  least  squares  (WLS)  framework.  The  WLS  framework 
is  general  and  subsumes  the  traditional  interpolation-based  theoretical 
foundation.  It  delivers  the  same  order  of  accuracy  but  far  superior  sta¬ 
bility  and  flexibility  on  irregular  meshes  over  complex  geometries,  a  fact 
which  we  have  proven  mathematically  and  demonstrated  numerically. 

Using  the  WLS  framework,  we  have  addressed  some  fundamental  prob¬ 
lems  in  scientihc  computing.  The  hrst  problem  is  surface  integration, 
which  is  a  core  procedure  for  a  variety  of  numerical  methods  such  as 
the  boundary  integral  method,  hnite  element  method,  surface  hnite  el¬ 
ements  etc.,  as  well  as  for  the  computation  of  geometrical  primitives 
such  as  surface  area  and  volume.  We  addressed  the  fundamental  ques¬ 
tion  of  whether  one  can  achieve  high-order  accuracy  given  a  piecewise 
linear  approximation  to  the  surface.  We  proposed  a  novel  method, 
which  can  compute  the  surface  integrals  to  high-order  of  accuracy  over 
piecewise-linear  meshes  [29,  9].  We  proved  theoretically  that  high-order 
of  accuracy  can  be  achieved  for  surface  integrals  using  a  least-  squares 
based  approximation  over  linear  meshes  and  demonstrated  numerical 
experiments  achieving  the  predicted  order  of  accuracy. 

2.2.  High-Order  Surface  Remeshing.  Surface  remeshing  is  a  widely 
used  strategy  for  a  variety  of  purposes  such  as  mesh  quality  improve¬ 
ment,  conformality  to  highly  curved  boundaries,  error-based  adapta¬ 
tion,  etc.  The  accuracy  of  geometry  can  be  easily  lost  if  low-order 
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point  projections  are  used  during  remeshing  operations.  We  investi¬ 
gated  the  problem  of  optimizing,  adapting,  and  untangling  a  surface 
triangulation  with  high-order  accuracy,  so  that  the  resulting  mesh  has 
sufficient  accuracy  to  support  high-order  numerical  methods,  such  as 
hnite  element  methods  with  quadratic  or  cubic  elements  or  generalized 
hnite  difference  methods.  The  fundamental  idea  to  solve  this  problem 
was  to  integrate  various  remeshing  strategies  with  a  high-order  point 
projection.  Specihcally,  we  developed  a  remeshing  strategy  based  on 
variational  mesh  optimization  along  with  high-order  surface  reconstruc¬ 
tion  that  produces  high-quality  triangular  meshes.  The  algorithm  was 
further  improved  to  allow  of  untangling  mildly  folded  triangles.  These 
results  were  reported  in  [4].  In  general,  high-order  point  projection 
works  well  for  sufficiently  resolved  regions.  However,  they  may  intro¬ 
duce  large  variations  if  applied  to  regions  which  he  beyond  the  asymp¬ 
totic  limit  of  polynomial  httings  such  as  under-resolved  regions  describ¬ 
ing  regions  of  high  curvature.  We  introduced  safeguards  (also  known  as 
limiters)  to  capture  such  problematic  regions  and  subsequently  reduce 
the  order  of  htting,  resulting  in  a  robust  algorithm  [28].  Extensive  nu¬ 
merical  experiments  were  performed  to  demonstrate  the  effectiveness 
of  our  method  in  obtaining  convergent  differential  quantities  as  well  as 
comparison  to  other  methods.  This  work  has  been  integrated  in  the 
front-tracking  package  FronTier. 

2.3.  Theoretical  Model  of  Curvature  Effect  of  Thin  Membranes 

As  part  of  the  project,  we  investigate  the  curvature  effect  of  a  thin, 
curved  elastic  interface  [47].  Such  an  interface  separates  two  subdo¬ 
mains,  such  as  the  canopy  of  parachutes,  and  it  exerts  a  pressure  due 
to  a  curvature  effect.  We  refer  to  this  pressure  as  the  interface  pressure; 
it  is  similar  to  the  surface  tension  in  fluid  mechanics.  It  is  important  in 
some  applications,  such  as  the  canopy  of  parachutes,  biological  mem¬ 
branes  of  cells,  balloons,  airbags,  etc.,  as  it  partially  balances  a  pressure 
jump  between  the  two  sides  of  an  interface.  In  this  project,  we  showed 
that  the  interface  pressure  is  equal  to  the  trace  of  the  matrix  product 
of  the  curvature  tensor  and  the  Cauchy  stress  tensor  in  the  tangent 
plane.  We  derive  the  theory  for  interfaces  in  both  2-D  and  3-D,  and 
present  numerical  discretizations  for  computing  this  quantity  over  tri¬ 
angulated  surfaces.  Figure  1  shows  an  example  computation  of  the 
curvature  effect,  where  a  torus  with  inner  and  outer  radii  0.3  and  1  m 
was  artihcially  expanded  into  one  with  inner  and  outer  radii  0.32  and 
1.05  m.  Colors  of  left  and  right  images  indicate  the  tangential  stress 
and  interface  pressure,  respectively.  It  can  be  seen  that  the  stress  (and 
also  pressure)  is  larger  at  the  outside  than  at  the  inside.  Figure  2  shows 
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the  high-order  convergence  rates  that  were  achieved  of  this  model  when 
discretized  using  the  weighted  least  squares  framework. 


stress 
mm  0.42 
0.41 
0.4 
0.39 


pressure 

H  0.95 
—  0.85 


Figure  1.  An  example  computation  of  interface  pres¬ 
sure  due  to  curvature  effect. 


Figure  2.  High-order  convergence  results  for  computed 
interface  pressure  under  grid  rehnement  for  the  torus. 


3.  Appligations 

3.1.  Brittle  Fracture.  We  have  developed  new  mathematial  models, 
numerical  algorithms  and  parallel  software  for  the  simulation  of  brittle 
fracture  of  solids.  The  algorithms  are  based  on  energy  minimization 
in  a  network  of  triangular  springs  with  critical  strain  and  splitting 
of  overstressed  bonds.  The  algorithms  ensure  conservation  of  mass 
during  the  crack  evolution.  Special  energy  penalty  terms  preclude  mesh 
folding  and  overlapping  in  the  simulation  of  compressed  materials.  The 
main  emphasis  of  the  research  was  on  the  study  of  brittle  fracture, 
but  elasto-plastic  models  for  springs  have  also  been  developed  for  the 
simulation  of  plastic  deformations  with  limited  shear  bands. 
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Figure  3.  Comparison  of  non-conservative  and  mass- 
conservative  methods  using  the  example  of  projectile  im¬ 
pact  induced  wall  fracture.  Left;  the  non-conservative 
method  results  in  large  mass  loss,  with  the  creatation  of 
a  large  void;  Right:  The  mass  conservative  method  leads 
to  increased  pressure  and  the  formation  of  cracks. 

Our  models  signihcantly  improved  similar  previous  models  based  on 
the  energy  minimization  of  spring  networks.  The  main  deficiences  of 
prior  methods,  resolved  by  our  methods,  are 

•  Algorithms  applicable  only  to  stretched  solids.  Unphysical  mesh 
folding  occurs  if  these  models  are  applied  to  solids  under  compres¬ 
sion  forces. 

•  Non-conservative  bond  removal  algorithm  which  leads  to  the  loss 
of  mass.  This  issue  is  less  important  if  only  a  small  number  of 
isolated  cracks  develop  in  streched  solids,  but  it  prevents  the  sim¬ 
ulation  of  comminuted  fracture  zones  in  compressed  solids. 

•  Inability  to  reproduce  complex  patterns  ocurring  in  brittle  frac¬ 
ture  processes  (for  instance  a  comminuted  zone  and  a  series  of 
radial  cracks  in  solid  plates  hit  by  high  velocity  projectiles). 

•  Only  2D  dimensional  models 

•  Serial  codes  for  running  on  single-processor  computers 

Our  representation  of  solids  by  spring  networks  contains  the  two  de¬ 
grees  of  freedom  necessary  to  match  real  material  properties  and  it 
exhibits  a  stable  Poisson  ratio.  Two  regimes  of  the  brittle  fracture 
have  been  considered:  adiabatically  slow  deformation  with  breakup. 
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Figure  4.  Brittle  fracture  of  thin  disks  hit  by  projec¬ 
tiles  exhibiting  comminuted  regions  and  radial  cracks. 
The  number  of  radial  cracks  and  the  comminuted  region 
size  decrease  with  the  increase  of  material  plasticity  and 
reduction  of  the  projectile  velocity. 


Figure  5.  Left  and  center:  2D  simulation  of  brittle  frac¬ 
ture  of  a  thick  wall  hit  by  a  projectile.  One-layer  of  in¬ 
scribed  circles  (left)  and  two-layers  of  inscribed  circles 
(center)  were  used  in  the  fragment  collision  detection  al¬ 
gorithm.  Right:  Fracture  of  a  3D  rod  under  tension. 

and  instantaneously  fast  deformation  with  the  formation  and  propaga¬ 
tion  of  cracks  in  stressed  materials.  Parallel  software  for  the  fracture 
of  brittle  materials  under  strain  has  been  developed  using  the  exter¬ 
nal  parallel  packages  TAO  and  Global  Arrays,  developed  within  DOE 
high  performance  computing  initiatives.  A  Schwartz-type  overlapping 
domain  decomposition  and  the  corresponding  acceleration  techniques 
have  also  been  studied.  Three  different  visualization  techniques  have 
been  developed  to  capture  details  of  fractured  zones  in  3D. 

The  software  has  been  applied  to  the  simulation  of  fracture  of  solids  in 
the  regime  of 
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1.  slow  stretching  deformations, 

2.  materials  posessing  high  initial  quasi-stable  overstress  such  as  in 
the  rapid  disintegration  of  highly  tempered  glass  in  the  phenom¬ 
enon  called  the  Prince  Rupert  drop,  and 

3.  fracture  of  thin  brittle  discs  and  thick  walls  hit  by  high  velocity 
projectiles  [48,  49]. 

The  bifurcation  of  the  fracture  dynamics  from  the  growth  of  the  com¬ 
minuted  zone  to  the  propagation  of  isolated  radial  cracks,  typical  for 
the  fracture  of  glass  sheets  and  thin  ceramic  plates  hit  by  projectiles, 
has  been  reproduced  in  our  numerical  experiments.  Scaling  studies  in¬ 
volving  the  variation  of  material  properties  and  projectile  velocity  have 
been  performed  (see  Figure  4).  2D  and  3D  simulations  of  the  fracture 
of  thick  walls  are  depicted  in  Figure  5  [48,  49].  The  fracture  model  has 
also  been  used  in  a  coupled  multiscale  simulation  of  the  nuclear  fuel 
rod  failure  within  a  study  of  nuclear  reactor  safety  issues  [2,  1]. 

The  main  research  goal  of  the  present  study  is  to  improve  the  avail¬ 
ability  of  the  method  by  use  of  a  hnite  elements  implemention.  This 
will  enable  standard  hnite  element  codes  to  simulate  complex  fracture 
regimes  and  hows  of  disintegrated  materials.  The  related  goal  is  the 
development  of  highly  scalable  software,  and  applications  to  problems 
relevant  to  Army  Research  Laboratory  and  DOE  Nuclear  Energy  re¬ 
search. 

3.2.  Numerical  Modeling  of  Parachute  Delivery  System.  Sim¬ 
ulation  of  parachute  inhation  via  computational  method  has  attracted 
the  attention  of  scientists  at  the  Department  of  Defense  laboratories 
and  academia  alike.  Some  of  the  successful  studies  include  those  re¬ 
ported  in  [34,  33,  30,  35,  31,  32],  where  the  Deforming  -  Spatial  -  Do¬ 
main  /  Stabilized  Space-Time  (DSD/SST)  method  [39,  40]  was  used 
in  [33,  30,  35,  31,  32].  The  studies  in  [43,  42,  37,  38,  44]  also  used 
the  DSD/SST  as  the  core  numerical  method,  but  involved  new  ver¬ 
sions  and  special  techniques.  These  studies  successfully  address  the 
computational  challenges  in  handling  the  geometric  complexities  of  the 
parachute  canopy  and  the  contact  between  parachutes  in  a  cluster.  Kim 
and  Peskin  et  ah  used  the  immersed  boundary  method  to  study  the 
semi-opened  parachute  in  both  two  and  three  dimensions  [18,  19],  their 
simulations  are  on  small  Reynolds  number  (about  300)  and  applied 
payload  with  several  grams.  Karagiozis  used  the  large-eddie  simula¬ 
tion  to  study  the  parachute  in  Mach  2  supersonic  flow  [16].  Purvis 
[26,  27]  used  springs  to  represent  the  structures  of  the  fore-body,  the 
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suspension  lines,  and  the  canopy,  etc.  In  these  papers,  the  authors 
used  cylindrical  coordinates  with  the  center  line  as  the  axis.  In  the 
paper  by  Strickland  et  ah  [36],  the  authors  developed  an  algorithm 
called  PURL  to  couple  the  structure  dynamics  (PRESTO)  and  fluid 
mechanics  (CURL),  in  which  mass  is  added  to  each  of  the  structure 
node  based  on  the  diagonally  added  mass  matrix  and  a  pseudo  is  com¬ 
puted  from  the  fluid  code  which  is  the  sum  of  the  actual  pressure  and 
the  pressure  associated  with  the  diagonally  added  mass.  Tutt  and  Tay¬ 
lor  [46,  45]  simulated  the  parachute  through  the  LS-DYNA  code.  They 
used  an  Eulerian-Lagrangian  penalty  coupling  algorithm  and  multi¬ 
material  ALE  capabilities  with  LS-DYNA  to  replicate  the  inflation  of 
small  round  canopies  in  a  water  tunnel. 

Our  computational  approach  to  the  simulation  of  the  parachute  deliv¬ 
ery  system  is  based  on  the  front  tracking  platform.  In  the  last  ARO 
funded  project  during  which  the  PI  served  as  the  principal  developer, 
we  established  a  computational  platform  for  the  parachute  study  by 
employing  the  spring  model  for  the  parachute  canopy  and  the  string 
chords.  We  designed  new  data  structure  to  allow  the  application  of  the 
FronTier  library  to  the  dynamic  motion  of  fabric  surface  driven  by  the 
gravitational  force  of  payload  and  the  fluid  pressure.  We  discretize  the 
fabric  surface  into  a  homogeneously  triangulated  surface  mesh. 


3.2.1.  Analysis  of  the  Spring  System.  We  made  a  detailed  analysis  of 
the  spring  system  which  is  used  to  model  the  fabric  surface  of  the 
parachute  canopy  (Figure  6).  When  no  external  driving  force  is  applied, 
the  fabric  surface,  which  is  represented  by  the  spring  mass  system,  is  a 
conservative  system  whose  total  energy  (kinetic  energy  plus  potential 
energy)  is  a  constant.  Assuming  each  mesh  point  represents  a  point 
mass  m  in  the  spring  system  with  position  vector  Xj,  the  kinetic  energy 
of  the  mass  point  i  is  where  Xj  is  the  time  derivative,  or 

velocity  vector  of  the  mass  point  i. 

We  considered  two  types  of  spring  systems.  The  first  one,  which  we 
refer  to  as  Model-I,  or  the  elastic  foil  model,  has  the  potential  energy 
between  two  mass  points  x*  and  Xj  in  the  form  of 

k 

^ij  =  2  -  xj)  -  (xio  -  Xjo)|^  ,  (1) 

where  k  is  the  spring  constant  and  Xjo  is  the  equilibrium  position  of 
mass  point  i.  This  potential  energy  leads  to  the  following  equation  of 
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Figure  6.  The  spring  model  on  a  triangulated  mesh. 
Each  vertex  point  in  the  mesh  represents  a  mass  point 
with  point  mass  m.  Each  edge  of  the  triangles  has 
an  equilibrium  length  set  during  initialization  and  the 
changing  length  exerts  a  spring  force  on  the  two  neigh¬ 
boring  vertices  in  opposite  directions.  Gores  are  added 
as  curves  with  larger  spring  constant.  This  plot  shows 
the  mesh  of  canopy  with  gores. 

motion  for  the  system 


where  x^o  is  the  equilibrium  position  of  mass  point  i.  This  system  can 
be  simplified  by  the  substitution  x'*  =  Xj  —  x^o,  i  =  1,2,  -  ■  ■  ,N,  which 
yield 


(3) 


This  system  does  not  match  the  properties  of  fabric,  but  it  can  be 
analyzed  in  closed  form. 

Using  the  Levy-Desplanques  Theorem  for  diagonally  dominated  ma¬ 
trix,  we  proved  [21]  that  the  motion  of  Model-I  is  purely  oscillatory 
and  there  exists  an  upper  bound  for  the  eigen  frequencies  of  Model-I. 
In  addition,  through  the  Gershgorin  circle  theorem,  we  have  obtained 


9 


an  upper  bound  for  the  eigen  frequency  of  the  oscillatory  motion 


u  < 


2Mk 


m 


(4) 


A  fabric  surface  is  considered  as  a  membrane  which  is  an  idealized  two 
dimensional  manifold  for  which  forces  needed  to  bend  it  are  negligible. 
Therefore  for  parachute  canopy,  we  should  not  consider  the  bending 
energy.  Model-I  contains  strong  bending  force  and  is  not  suitable  for 
fabric  modeling.  To  obtain  a  realistic  spring  system  for  the  simulation 
of  fabric  surface,  we  have  to  assume  that  the  spring  force  between 
two  neighboring  point  masses  is  only  proportional  to  the  displacement 
from  their  equilibrium  distance,  the  potential  energy  due  to  the  relative 
displacement  between  two  neighboring  point  masses  x*  and  Xj  is  given 
by 

Idj  =  ^(|xi-xj| (5) 

where  =  |xjo  —  x^qI  is  the  equilibrium  length  between  the  two  point 
masses.  Here  we  have  made  a  modihcation  to  the  model  used  by  Choi 
and  Ko  in  that  we  allow  both  contraction  and  expansion  forces  while 
in  Choi  and  Ko’s  equation,  no  force  is  applied  if  |xj  —  x,|  <  1%.  Choi 
and  Ko’s  choice  does  not  conserve  energy  and  would  allow  a  surface  to 
shrink  to  a  point  without  resistance.  Such  shrinking  is  unrealistic  for 
a  fabric  surface  such  as  the  parachute  canopy.  We  refer  to  this  model 
as  Model-II,  or  the  fabric  model. 


VVllPP^ll  cxo 


The  equation  of  motion  of  Model-II  can  be 

^  =  -Z]«u(Cj)(xi-x,), 


m 


(6) 


where 

aij  =  Vijk 

V  '  q 

Upon  linearization,  we  have 

Using  Bij  to  make  dot  product  for  each  term,  we  obtain  the  tangential 
component  of  acceleration 

fk  =  -Bij  ■  riijk  ®  Gij  +  ~ 

=  -rjijk  {Bij  ■  ((5xi  -  5xj)) 


(8) 
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This  equation  indicates  that  the  tangential  force  of  the  fabric  model 
(Model-II)  is  the  same  as  in  Model-I.  Therefore  we  conclude  that  for 
Model-II,  only  the  tangential  motion  is  oscillatory. 


Figure  7.  Convergence  test  of  the  drum  membrane  un¬ 
der  mesh  rehnement.  From  left  to  right  the  computa¬ 
tional  mesh  of  the  domain  are  15,  30,  and  60  respectively. 
The  total  mass  of  the  membrane  is  kept  constant  in  the 
simulations.  The  upper  three  plots  show  the  membrane 
position  at  t  =  1  and  the  lower  plots  show  the  membrane 
position  at  t  =  2. 


3.2.2.  Numerical  Convergence  Verification.  Numerical  convergence  un¬ 
der  computational  mesh  rehnement  is  a  crucial  step  in  assessing  the 
validity  of  spring  mass  model.  For  this  purpose,  we  carried  out  two 
sequences  of  simulations  with  increasing  number  of  grid  points.  Our 
convergence  study  includes  both  one  dimensional  string  in  two  dimen¬ 
sional  plane  and  two  dimensional  membrane  in  three  dimensional  space. 

In  one  case  of  the  verihcation  study,  a  circular  vibrating  membrane  with 
radius  of  0.4m  and  total  weight  of  dSOgf  is  considered.  The  membrane  is 
linearly  perturbed  at  initial  time  and  its  circular  boundary  is  hxed.  A 
sequence  of  four  cases  with  computational  mesh  15^,  30^, 60^, and  120^ 
are  computed.  The  errors  of  total  area  are  shown  in  Table  1.  From 
Table  1,  we  can  clearly  see  that  the  Cauchy  errors  are  decreasing  as  the 
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reference  mesh  size 

ca 

Cp 

15  and  30 

0.02528 

1.91810 

1.97967 

30  and  60 

0.01507 

1.21784 

1.21772 

60  and  120 

0.00604 

0.53550 

0.52837 

Table  1.  Convergence  tests  of  spring  model  for  a  fabric 
drum.  In  the  computational  sequences,  the  total  mass 
of  the  membrane  is  hxed.  As  the  number  of  points  in¬ 
creases,  the  point  mass  is  reduced  accordingly.  Cauchy 
error  is  calculated  on  two  consecutive  mesh  sequences. 
Column  ca,  Cfc,  and  Cp  are  errors  of  total  area,  total  ki¬ 
netic  energy  and  total  spring  potential  energy  respec¬ 
tively.  The  numerical  results  show  the  hrst  order  conver¬ 
gence. 


computational  mesh  is  rehned.  The  convergence  rate  for  the  membrane 
is  also  of  hrst  order. 


3.2.3.  Stress  Analysis.  Stress  analysis  on  the  parachute  canopy  and 
the  string  chords  during  the  process  of  inhation  is  very  important  to 
the  held  test  of  the  parachute.  The  natural  stress  on  each  side  of  the 
triangle  in  the  spring  model  is  the  restoring  force  due  to  the  stretching 
of  each  triangle  side.  Let  Ti,r2,r3  be  the  natural  stress  on  side  1,2,3 
of  a  triangle  respectively,  we  have 


T,  =  k(l,  -  If),  i  =  1,2,3, 


(9) 


where  k  is  the  spring  constant,  is  the  equilibrium  length  of  side  i 
and  li  is  the  stretched  length  of  the  side  i.  This  natural  stress  can  be 
converted  into  the  stress  in  Cartesian  coordinates  on  the  plane  of  the 
triangle.  The  Cartesian  stress  is  a  2  x  2  tensor  in  the  plane  of  the 
triangle 


a  =  f 

The  conversion  from  natural  stress 
mapping  matrix,  that  is 


(^xy  I 

'^yy  ) 

to  Cartesian 


(10) 


stress  is  through  a 


SiCi 

-S2C2 

-S3C3 


(11) 


where  c*  =  cos6*j  =  dxi/li  and  Si  =  sin6*j  =  dyi/k  are  the  trig  func¬ 
tions  of  the  angle  between  each  side  and  the  x— axis.  The  stresses  in 
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two  principal  directions  are  obtained  via  diagonalization  of  the  stress 
tensor,  that  is 


^XX 

^xy 


xy 

yy 


=  T 


-1 


CTl 

0 


0 

0-2 


T, 


(12) 


where  cxi  and  a2  are  the  solutions  of  the  characteristic  equation 


^xx  ^1,2  ^xy 

^  xy  ^yy  ^1,2 


0. 


(13) 


We  use  the  von  Mises  stress 

^vm  =  yaf  +  Cr|  -  CTiO-s  (14) 

to  measure  the  tension  on  the  fabric  material  surface.  The  safety  fac¬ 
tor  of  the  material  is  defined  as  the  ratio  of  the  significant  strength  of 
the  material  to  the  von  Mises  stress.  When  this  factor  is  decreased  to 
a  critical  value  (which  corresponding  to  high  stress  of  the  fabric  sur¬ 
face),  it  sends  a  warning  signal  for  the  design  of  the  parachute  canopy. 
Figure  8  shows  the  von  Mises  stress  in  color  in  three  fabric  surface 
simulations. 


Figure  8.  The  von  Mises  Eq.  (14)  stress  in  the  spring 
model.  The  left  plot  shows  the  von  Mises  stress  of  a 
rectangular  membrane  when  pulled  from  the  four  cor¬ 
ners.  The  center  plot  shows  the  von  Mises  stress  of  the 
triangular  membrane  pulled  from  the  three  vertices.  The 
right  plot  is  the  stress  on  parachute  canopy  after  infla¬ 
tion. 


3.2.4.  Fluid  Canopy  Coupling.  To  model  the  parachute  system,  an  ac¬ 
curate  coupling  between  the  Naiver-Stokes  equation  and  the  structure 
dynamics  must  be  carefully  considered  near  the  canopy  surface.  The 
method  we  designed  for  the  simulations  of  air  delivery  system  uses  the 
superposition  of  impulse  on  every  mass  point.  Each  mass  point  in  the 
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spring  system  acts  as  an  elastic  boundary  point  and  exerts  an  impulse 
on  the  fluid  in  its  normal  direction.  Our  algorithm  ensures  that  the 
action  and  reaction  between  the  spring  mass  point  and  the  fluid  solver 
are  equal  in  magnitude  and  opposite  in  directions,  a  requirement  of 
Newton’s  third  law. 

We  applied  Peskin’s  delta  function  immersed  boundary  method  [18,  19] 

f(x,  t)  =  J  F{s,t)5{x  — 'K.{s,t))ds.  (15) 

The  difference  between  our  method  and  Kim  and  Peskin’s  method  lies 
in  the  calculation  of  F.  Instead  of  computing  the  tension  through  the 
derivative  with  respect  to  the  arc  length,  we  use  the  impulse  of  the 
mass  point  as  a  result  of  the  superposition  of  three  forces  from  the 
spring  system,  that  is 

F(xj,i)  =  d(Ig  +  Ip  +  IJ/dt  (16) 

Eq.  (16)  is  more  physically  realistic,  especially  because  Is  is  obtained 
from  the  spring  equations.  In  the  canopy  spring  system,  we  have  ob¬ 
served  that  the  tension  is  high  at  the  top  of  the  canopy  where  the 
curvature  is  almost  zero. 


3.2.5.  Validation  of  the  System.  Our  primary  objective  is  to  conduct 
the  fabric-fluid  coupled  simulations  and  validate  it  against  the  available 
parachute  experiments.  Parachutes  differ  in  geometry  and  dimensions 
of  the  canopies  and  risers.  Among  these  parachutes,  the  T-10  parachute 
is  used  by  the  Army  as  a  personnel  carrier.  This  type  of  parachute  has 
a  parabolic  shape  for  its  canopy  with  a  vent  at  the  top.  The  G-11 
parachute  is  a  cargo  parachute  with  no  vent.  Its  dimensions  could  vary 
and  is  usually  used  as  a  multiple  parachute  system  to  deliver  support¬ 
ing  equipment.  We  have  studied  the  single  1/3-scale  G-11  parachute 
inflation.  The  cross  parachute  has  four  open  side  vents  and  can  be 
used  as  a  sports  parachute  and  small  cargo  carrier.  Figure  9  shows  the 
inflation  sequence  of  a  cross  parachute.  These  three  parachutes  pro¬ 
duce  different  patterns  of  airflow  around  the  canopy  and  exert  different 
drags  to  the  parachute  system. 

For  comparison  with  the  indoor  vertical  parachute  test,  we  initialized 
the  parachute  with  flat  circular  canopy  of  2.134  m  (7  ft)  nominal  diam¬ 
eter.  We  carried  out  pre-step  running  to  deform  the  parachute  canopy 
and  then  used  the  resulting  geometry  as  the  initial  state  for  the  drop 
test  simulation.  The  initial  skirt  diameter  makes  different  time  graph 
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for  full  inflation.  In  this  test,  there  is  a  breathing  motion  or  the  over¬ 
inflation  of  the  canopy.  It  is  shown  that  the  canopy  progresses  into  an 
over-inflated  shape,  almost  flat,  and  then  contracts  and  over-inflates 
again.  The  breathing  motion  is  an  oscillatory  motion.  The  canopy 
appeared  to  expel  the  excess  air  by  means  of  the  breathing.  In  the 
experiment,  the  breathing  motion  was  also  caused  by  the  constraint  on 
the  parachute,  imposed  by  the  guide  wire.  The  breathing  motion  in  the 
simulation  is  smaller  because  there  is  no  vertical  motion  restriction  such 
as  a  guide  wire.  The  terminal  velocity  of  canopy  has  a  good  agreement 
with  the  vertical  parachute  test  results.  The  experimental  data  shows 
that  the  descent  speed  rises  rapidly  to  a  peak.  It  slows  down  while  the 
parachute  inflates  and  then  slowly  approaches  a  steady  descent  speed 
of  4.27  mjs  (14  ft/s).  Even  though  the  numerical  solution  cannot  be 
compared  with  the  experimental  data  during  the  initial  period  of  time 
in  the  simulation,  we  have  observed  that  it  reaches  the  terminal  ve¬ 
locity  of  about  3.9  m/s.  The  difference  between  the  terminal  velocity 
in  numerical  simulation  and  the  experimental  data  is  thought  due  to 
the  lack  of  porosity  of  the  canopy  in  the  current  numerical  model.  We 
also  simulated  the  1/3-scale  G-11  parachute  with  different  canopy  fold¬ 
ing  level,  they  all  converge  to  the  terminal  speed  of  3  m/s  in  several 
seconds. 

Parachute  breathing  is  an  oscillatory  motion  caused  by  the  interaction 
between  the  fluid  force  and  the  canopy  at  the  skirt  of  the  parachute. 
The  large  difference  between  upper  and  lower  side  pressure  causes  vi¬ 
bration  in  the  projected  drag  area  of  the  canopy  and  thus  the  oscillation 
in  descent  speed.  Numerical  study  of  breathing  frequency  is  important 
for  understanding  the  stability  of  the  parachute.  In  our  simulation, 
we  initialized  2.134  m  (7  ft)  flat-circular  parachute  to  compare  with 
experimental  data  by  Tutt  et  ah  [45].  The  breathing  frequency  in  our 
simulation  is  approximately  1.6  -  2.0  Hz  (0.5  s  to  0.6  s  period).  This 
is  in  the  acceptable  range  with  the  experimental  frequency  which  is 
2.0  Hz  (0.5  s  period).  If  we  use  a  parachute  with  larger  size  or  lower 
Reynolds  number,  the  breathing  period  will  increase.  For  example,  the 
average  period  of  breathing  for  the  10.7  m  T-10  parachute  is  2.3  s  in 
experiment  [10],  while  it  is  about  2  s  in  our  simulation. 


3.2.6.  New  Computational  Development.  In  the  last  few  months,  we 
have  accomplished  several  new  tasks  in  the  parachute  simulation  code 
which  enabled  us  for  more  realistic  and  robust  computation  on  the 
parachute  problem. 
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Figure  9.  Simulation  of  cross  parachute  unfolding  and 
inflation.  The  starting  state  of  the  parachute  is  partially 
folded  state.  The  parachute  has  a  small  horizontal  drift 
because  the  folded  state  is  not  perfect  symmetric. 


Parachute  gores  are  stitched  by  the  reinforcement  cables.  The  re¬ 
inforcement  cables  are  important  structures  on  the  parachute  canopy 
surface.  To  accurately  model  the  aerodynamic  motion  of  a  parachute, 
we  need  to  create  a  mathematical  model  which  reveals  the  geometry 
as  well  as  the  material  strength  of  the  canopy  surface  and  the  gores. 
The  gore  structures  in  the  parachute  system  have  important  effect  on 
the  stability  of  the  parachute  motion.  In  our  model,  we  treat  the  rein¬ 
forcement  cables  as  curves  embedded  in  the  canopy  surface.  Young’s 
modulus  for  the  fabric  surface  and  the  reinforcement  cable  are  set  at 
different  values.  This  is  realized  by  assigning  different  spring  constants 
to  the  surface  mesh  and  to  the  gore  curves.  The  insertion  of  gores  as 
stiffened  interior  curves  in  the  canopy  surface  mesh  is  demonstrated  by 
Figure  6.  Figure  10  shows  that  when  fully  inflated,  the  spring  system 
reveals  both  the  patches  and  the  indentation  of  the  reinforcement  ca¬ 
bles.  The  effect  of  the  gores  on  suppressing  the  vortex  flow  at  the  top 
of  the  canopy  is  being  studied.  This  will  be  published  in  a  new  paper 
we  are  currently  working  on. 

Fabric  porosity  is  another  important  property  which  must  be  cor¬ 
rectly  modeled  for  the  canopy  surface.  It  has  long  been  known  that 
fabric  permeability  is  an  important  factor  in  parachute  design.  By 
carefully  considering  the  balance  between  payload  and  drag,  a  hnite 
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Figure  10.  The  reinforcement  cables,  or  the  gores,  are 
modeled  using  the  same  spring  system,  but  with  differ¬ 
ent  spring  constant  along  the  gore  curves.  The  left  plot 
shows  fully  opened  canopy  with  gores.  After  the  infla¬ 
tion,  the  gore  structure  is  clearly  revealed.  The  right  plot 
shows  the  surface  mesh  and  the  details  of  the  gore  curves. 


Figure  11.  Angled  deployment  of  C-9  parachute  with 
the  flow.  The  deployment  angle  between  the  initial 
parachute  and  the  direction  of  flow  are  (from  left  to  right) 
15°,  30°,  45°,  60°,  75°,  respectively. 


permeability  will  make  the  parachute  much  more  stable.  The  perme¬ 
ability  of  the  parachute  material  often  plays  a  vital  role  in  the  design 
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of  parachute.  A  state-of-art  tuning  between  an  impervious  material 
to  a  highly  permeable  fabric  can  make  a  parachute  from  a  wandering 
sloth  into  a  plummeting  stabilizer.  We  have  considered  three  different 
implementations  of  porosity  for  the  parachute  including  Tezduyar  and 
Sathe  [41],  Kim  and  Peskin  [18],  and  Tntt[46].  We  have  chosen  the 
latter  two  implementations  because  they  use  the  similar  frame  work 
of  computational  mesh  as  we  do  and  are  easy  to  adapt  into  onr  front 
tracking  environment. 


For  Kim  and  Peskin’s  method,  porosity  is  considered  as  the  leaking 
pores  in  the  immersed  elastic  bonndary.  Let  (3  be  the  density  of  pores, 
which  means  there  are  j3ds  pores  in  the  interval  (s,  s  +  ds)  along  the 
snrface.  If  each  pore  has  an  aerodynamic  conductance  eqnal  to  7,  which 
means  the  flux  through  the  pore  is  7(^1  —  P2)  where  pi  and  p2  are  the 
pressnres  on  the  two  sides  of  the  boundary,  then  the  flux  through  the 
interval  (s,  s  +  ds)  of  the  bonndary  is  given  by  P'yipi  —  P2)ds. 


Tntt’s  algorithm  is  closest  to  the  front  tracking  implementation.  Tutt 
used  the  Ergun  equation  to  describe  the  magnitude  of  porous  flow 
velocity  at  a  given  pressure  difference  based  on  two  coefficients  in  the 
following  eqnation: 


AP 

~T~ 


P  2  T  2 


where: 


K, 


150  •  (1  -e)2 


^  1.75 -(l-e) 

are  referred  to  as  the  viscous  and  inertial  factors  respectively.  D  is 
dehned  as  the  characteristic  length,  e  is  the  porosity  and  is  equal  to 
the  ratio  of  the  void  and  total  volume.  Vf,  p  and  p  are  flnid  velocity, 
viscosity  and  density  respectively.  In  our  FronTier-airfoil  code,  the 
pressure  difference  is  compnted  after  projecting  the  velocity  into  the 
divergence-free  space.  We  can  then  use  the  pressure  difference  on  two 
sides  of  the  canopy  to  compute  the  Vf  and  add  it  to  the  advection  solver 
as  flux  from  the  immersed  elastic  surface. 


Collision  handling  on  parachnte  contact  is  throngh  the  modified  Fron 
Tier  library.  The  computational  module  for  parachute  simulation  is 
built  on  the  FronTier  library.  We  have  used  many  existing  data  strnc- 
tures  and  functionality  for  geometrical  handling.  However,  the  study  of 
parachute,  airfoil  and  other  fluid  structure  based  simulations  reqnires 
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major  revisions  to  the  software  library.  Among  those  needed  are  fnnc- 
tions  to  handle  the  non-manifold  snrface  and  three  dimensional  curves 
(not  as  boundary  of  a  surface,  such  as  the  string  chords).  The  original 
FronTier  code  was  designed  for  the  study  of  fluid  interface  instabil¬ 
ity  problems.  For  these  problems,  collision  surfaces  are  resolved  by 
merging  or  bifurcation.  However,  fabric  surface  can  neither  merge  nor 
bifurcate,  therefore  we  need  to  have  functions  which  can  carefully  deal 
with  the  repulsive  contact  and  collision. 

Global  indexing  is  a  new  feature  added  to  the  parachute  code.  The 
original  FronTier  had  to  deal  with  frequent  surface  mesh  optimiza¬ 
tion  and  topological  reconstructions.  This  makes  the  parallelization 
based  on  global  index  very  difficult.  As  a  result,  the  original  FronT 
ier  library  relied  on  floating  point  matching  for  parallel  communica¬ 
tion.  The  floating  point  matching  is  not  100  percent  reliable  and  some 
more  complicated  algorithms  have  to  be  implemented  as  reinforcement. 
However  for  fabric  surface,  especially  when  spring  model  is  used,  the 
inter-connectivity  and  proximity  of  the  interface  marker  points  are  not 
changed.  Therefore,  global  indexing  is  ideal  for  the  parallel  commu¬ 
nication  of  the  surface.  We  have  added  a  new  parallel  code  for  this 
functionality. 

Modularization  is  emphasized  in  our  code  development.  The  parachute 
module  is  an  independent  application  program.  This  new  module  con¬ 
sists  of  four  components:  (1)  the  initialization  module,  (2)  the  ODE 
module  for  the  spring  system,  (3)  the  PDE  module  for  fluid  dynamics, 
and  (4)  the  FronTier  library  for  the  interface  geometry  handling.  A 
total  of  about  15,000  lines  of  code  have  been  written  for  the  first  three 
component.  An  additional  4,500  lines  have  been  added  to  FronTier  to 
adapt  the  library  for  the  computation  of  fabric  surface. 

GPU  technology  [20]  is  to  use  the  Graphics  Processing  Unit  (GPU) 
together  with  the  GPU  to  accelerate  a  general-purpose  scientihc  and  en¬ 
gineering  application.  GPU  computing  can  offer  dramatically  enhanced 
application  performance  by  offloading  computation-intensive  portions 
of  the  programming  code  to  the  GPU  units,  while  the  remainder  of  the 
code  still  runs  on  the  GPU. 

Joint  GPU/GPU  application  is  a  powerful  combination  because  GPUs 
consist  of  a  few  cores  optimized  for  serial  processing,  while  GPUs  con¬ 
sist  of  thousands  of  smaller,  more  efficient  cores  designed  for  massive 
parallel  calculations.  Serial  portions  of  the  code  with  intensive  logi¬ 
cal  comparisons  run  on  the  GPU  while  floating  point  intensive  parallel 
portions  of  the  code  run  on  the  GPU. 


19 


3.2.7.  Delingette  Modification.  Recently,  Delingette  [5]  has  studied  spring 
model  on  triangulated  mesh  on  nonlinear  membrane.  This  new  model 


is  favored  because  of  its  derivation  from  the  elastic  membrane  model 


in  continuum  mechanics. 


Figure  12.  (a)  Rest  Triangle  T^o  whose  vertices  are 

(b)  Deformed  triangle  Tx  whose  vertices  are  x* 

The  idea  of  Delingette’s  revised  model  can  be  demonstrated  by  Fig¬ 
ure  12.  The  energy  of  membrane  W{T^o)  required  to  deform  a  single 
triangle  T^o  with  vertices  {x^,  x^,  Xg}  into  its  deformed  position  Tx  with 
vertices  {xi^,X2,X3}  consists  of  two  parts, 

•  The  energy  of  three  tensile  springs  that  prevent  edges  from  stretch¬ 
ing. 

•  The  energy  of  three  angular  springs  that  prevent  any  change  of 
vertex  angles. 

We  use  A^o,  ai  and  to  represent  the  initial  area,  initial  angles  and 
initial  lengths  of  the  rest  triangle  TxO.  Correspondingly,  we  use  ^x,  A) 
and  lij  to  denote  the  area,  angles,  and  lengths  of  the  deformed  triangle 
Tx.  The  edge  elongation  can  be  written  as  dUj  =  kj  —  Based  on  the 
analysis  by  Delingette  [5],  we  have 


(/°fc)^(2cot^aj(A  +  ^)  +  li) 


SAxO 


where 
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is  the  tensile  stiffness  and 


cot  ai  cot  aj{\  +  /x)  -  /i) 

8^1x0 


is  the  angnlar  stiffness.  7  and  /i  are  the  lame  coefficients  of  the  material. 
Those  coefficients  are  simply  related  to  the  two  physically  meaningfnl 
parameters  dehned  in  planar  elasticity  for  a  membrane,  that  is,  the 
Yonng  modulus  E  and  the  Poisson  coefficient  v  [6]: 

Ev  ,  Eil  —  V 
A  =  - - ^  and  /i  =  — - ^ 


The  Young  modulus  quantihes  the  stiffness  of  the  material,  whereas 
the  Poisson  coefficient  characterizes  the  material  compressibility. 


Applying  the  Rayleigh-Ritz  analysis,  it  is  considered  that  the  fabric 
surface  represented  by  triangular  mesh  should  evolve  by  minimizing 
its  membrane  energy,  therefore  along  the  opposite  derivative  of  that 
energy  with  respect  to  the  nodes  of  the  system,  that  is,  the  deformed 
positions  xp 


F^{T^o)  = 


dW{T^ 

(9x,' 


Ro 


{dl 


ik ) 


^ik 


Hk 


The  membrane  energy  to  deform  the  whole  triangulation  is  simply  the 
sum  of  the  energy  of  each  triangle.  Thus,  if  an  edge  is  shared  by  two 
triangles  Ti  and  T2,  the  total  stiffness  of  that  edge  k'^  is  the  sum  of  the 
two  stiffness  coefficients:  k"’"  =  k'^^  + 


If  the  second  term  in  Eq.  (17)  is  neglected,  Delingette’s  model  is  the 
same  as  our  original  model  except  that  the  spring  constant  will  vary  if 
the  corresponding  initial  triangle  deviate  from  the  isosceles  triangle. 


3.3.  Turbulence  Modeling. 


3.3.1.  Theory.  We  exploited  systematically  the  K41  theory  of  turbu¬ 
lence.  This  description  of  the  velocity  statistics  is  in  effect  a  Sobolev 
inequality,  and  assuming  this,  we  derive  an  existence  theorem,  not  just 
for  the  Navier-Stokes  equation,  but  for  the  limiting  inhnite  Reynolds 
number  Euler  equation  [3].  Pursuing  this  same  set  of  scaling  ideas, 
we  interpreted  the  derivation  of  subgrid  terms  in  terms  of  the  renor- 
maliation  group,  with  higher  order  corrections  [7].  Perhaps  the  most 
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important  idea  to  emerge  from  this  line  of  reasonaing  is  a  clear  pic¬ 
ture  of  physically  based  nonuniqueness  of  solutions  for  the  LES  tur¬ 
bulence  and  turbulence  mixing  simulations  and  for  the  limiting  Euler 
equations‘ (whose  solution  we  just  established).  The  idea  is  that  the 
compressible  Euler  equation  and  the  multifluid  Euler  equation  (com¬ 
pressible  or  not)  have  dimensionless  parameters  related  to  the  Schmidt 
and  Prandtl  numbers  and  to  the  ratio  of  the  bulk  to  shear  viscosity. 
For  LES  simulations,  the  parameters  have  turbulent  versions,  not  spec¬ 
ified  by  laws  of  physics,  but  rather  by  the  requirement  that  the  coarse 
grid  simulations  should  converge  as  best  as  possible  to  highly  resolved 
(inhnite  Reynolds  number)  simulations.  Since  the  latter  can  not  be 
achieved  in  practice,  the  requirement  is  in  effect  unspecified,  and  is 
only  settled  by  comparison  to  experiment.  These  parameters  intro¬ 
duce  nonuniqueness  to  the  LES  solutions.  This  is  in  some  sense  a  crisis 
for  applied  computational  science,  which  generally  believes  that  a  nu¬ 
merically  correct  simulation  will  give  rise  to  a  physically  correct  answer. 
We  see  that  this  statement  is  only  asymptotically  correct  in  the  limit 
of  inhnite  mesh  resolution,  that  is  for  direct  numerical  simulations.  To 
repair  this  crisis,  we  introduced  a  validation  program,  which  was  quite 
successful,  relative  to  the  range  of  experimental  data.  The  huid  mixing 
data  set  runs  up  to  Re  =  35,000,  but  many  important  problems  arise 
at  higher  Reynolds  numbers.  For  this  regime,  we  demonstrated  that 
a  numerical  extrapolation  from  the  experimental  Reynolds  numbers  to 
much  higher  ones  is  quite  mild,  so  that  effective  validation,  with  a  mild 
layer  of  extrapolation,  i.e.,  verihcatoin,  will  suffice  [24,  25]. 

In  current  and  ongoing  work  we  have  begun  the  extension  of  these 
results  to  the  convergence  of  the  cumulative  distribution  functions  (the 
hrst  integral  of  the  probability  density  functions).  A  statistical  tool, 
called  W*,  was  developed  to  assist  in  the  statistical  data  analysis,  for 
this  modihed  notion  of  stochastic  convergence  [17]. 

The  stocahstic  convergence  was  demonstratred  in  a  series  of  papers 
[14,  15,  22],  and  continues  with  our  current  research.  See  Figs.  13  and 
14. 


3.3.2.  Simulations.  [8,  11]  [12]  [13]  The  observation  of  Sec.  3.3.1  that 
turbulent  mixing  simulations  and  all  but  the  simplest  (incompressible, 
constant  density  single  fluid)  LES  turbulent  simulations  are  physically 
underspecihed  and  possess  nonunique  solutions  has  profound  implica¬ 
tions  for  numerical  simulations.  Even  this  one  exceptional  case  has 
been  shown  to  have  nonunique  solutions  in  a  mathematical  analysis. 
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Figure  13.  Spieces  concentration  in  a  flnid  mixtnre 
taken  in  a  slice  throngh  the  middle  of  the  mixing  zone. 
Left  to  right:  coarse,  mednin  and  hne  grids.  Visnal  in¬ 
spection  snggests  some  kind  of  mesh  convergence  of  the 
mixtnre  fractons,  from  a  statistical  point  of  view.  See 
Fig.  14  for  a  qnnatitative  measnre  of  this  convergence. 


Becanse  of  the  nonnniqneness,  verification,  that  is  convergence  nnder 
mesh  rehnement,  is  insnfficient  to  assure  physically  correct  answers. 
Validation,  that  is  comparison  to  experiment,  is  essential  for  scientihc 
accuracy.  Fortunately,  the  algorithm  we  developed.  Front  Tracking 
combined  with  subgrid  scale  terms  (FT/SGS)  gives  basically  exact  re¬ 
sults  compared  to  experiment  for  LES  turbulent  mixing  [8,  11].  See 
Fig.  16  for  an  example  of  solution  validation.  These  two  ingredients. 


23 


0.00  0.30  O.SO  0,00  0.30  0 ,  SO 

Figure  14.  Li  norm  differences  of  the  cumulative  dis¬ 
tribution  functions,  from  the  concentrations  showm  in 
Fif.  13.  The  two  plots  compare  coarse  to  hne  (left)  and 
medium  to  hne  (right).  The  decrease  error  measured  in 
the  Li  norm  under  mesh  refienment  is  a  quantitiative 
measure  of  the  stochastic  convergence. 
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Figure  15.  LES  simulations  for  a  shock  wave  induced 
mixing  problem  using  two  different  choices  of  the  sub¬ 
grid  model  coefficients.  The  difference  in  the  solutions 
indicates  the  nonuniqueness  of  the  simulation. 


tracking  and  subgrid  models,  play  different  and  opposite  roles.  The 
SGS  terms  add  turbulent  diffusion  and  viscosity,  while  the  front  track¬ 
ing  prevents  excess  thermal  and  species  diffusion.  We  believe  that  any 
good  numerical  scheme  with  both  of  these  elements  will  obtain  a  nearly 
correct  answer.  However,  final  tuning  to  account  for  algorithm  depen¬ 
dent  changes  still  leads  to  required  adjustments  to  achieve  validation 
goals.  The  tuning  will  amount  to  adjustment  of  diffusion  terms  [24]. 

In  Fig.  15,  we  show  the  probability  density  functions  from  two  simula¬ 
tions  of  a  shock  driven  turbulent  mixing  simulation,  but  with  different 
subgrid  model  coefficients.  The  disparity  of  the  solutions  is  a  numerical 
confirmation  of  nonuniqueness. 

Uncertainty  quantification  was  addressed  in  these  studies.  Our  valida- 
ton  studies  concern  the  edge  of  the  mixing  zone  for  Rayleigh-Taylor 
instabilities.  Unmeasured  long  wave  lengths  in  the  initial  conditions 
were  assessed  from  early  time  experimental  plates,  and  reconstructed, 
with  quantified  uncertainty,  leading  to  an  overall  uncertainty  in  the 
hnal  experimental  value  of  ±5%  [8,  11],  see  Fig.  16. 

Current  and  future  work  will  focus  on  the  convergence  of  the  probabil¬ 
ity  distribution  functions  for  the  mixtures.  Initial  results  show  evidence 
of  such  convergence  for  the  first  integral  of  the  pdfs,  namely  the  cu¬ 
mulative  distribution  functions  [25,  13];  These  same  references  display 
evidence  for  the  nonuniqueness  of  the  solutions  as  the  SGS  terms  (or 
equivalently  the  algorithm)  are  varied. 
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Figure  16.  A  comparison  study  of  three  simulations 
and  experimental  data  points,  for  an  acceleration  driven 
fluid  instability  problem.  The  three  simulations  have  the 
nominal  long  wave  length  (reconstructed  from  experi¬ 
ment)  initial  conditions,  and  the  same  multiplied  by  0 
or  2,  with  this  (large)  margin  set  from  uncertainties  in 
the  methods  of  data  analysis.  The  observed  variation  of 
the  instability  growth  rate  is  ±5%. 

3.3.3.  Tech  Transfer.  LANL  We  simulate  turbulent  mixing  in  inertial 
conhnement  flows.  As  a  base  calculation,  we  hnd  simulation  models 
which  show  improved  ht  to  experimental  data.  These  models  will  be 
used  to  test  the  turbulent  mixing  and  stochastic  convergence  models 
developed  here,  in  an  A  B  comparison  (with  and  without  enhanced 
resolution  methods. 

ORNL  We  study  a  highly  agitated  two  phase  flow  designed  to  enhance 
chemical  reactions  occuring  at  the  interface  between  the  two  phases 
[23].  A  correction  to  the  Navier-Stokes  equation,  an  added  pressure. 


26 


called  disjuncture  pressure  a  surface  effect  similar  in  origin  to  surface 
tension))  is  needed  to  capture  experimentally  observed  small  phenom¬ 
ena  of  local  phase  separation.  The  extra  pressure  is  a  surface  related 
force  familiar  in  lubrication  theory,  and  has  an  origin  in  thermodynam¬ 
ics  similar  to  that  of  surface  tension.  It  retards  the  merging  of  bubbles 
embedded  in  the  two  phase  flow  and  greatly  enhances  the  surface  area, 
thereby  accelerating  the  chemical  reactions. 

Stanford  University  We  study  turbuent  combustion  in  the  engine  of 
a  scram  jet.  The  main  thrust  of  the  effort  was  to  develop  hnite  rate 
chemistry  simulation  models.  We  also  study  interaction  of  turbulence, 
particle  laden  flow  and  radiation.  Depending  on  the  Stokes  number, 
the  particles  cluster  in  narrow  regions  of  high  strain  rate,  exterior  to 
vortex  tubes. 

4.  Impact  to  Education  and  Training 

The  previous  research  grant  also  supported  Joungdong  Kim,  a  perma¬ 
nent  resident  of  the  United  States.  Parachute  simulation  is  the  main 
topic  of  Dr.  Kim’s  PhD  thesis.  After  graduation,  Joungdong  Kim 
joined  the  Mathematics  Department  of  University  of  Texas  A&M  as  a 
visiting  assistant  professor.  Another  former  student  from  our  research 
group,  Eric  Rathnayke  joined  Military  Accessions  Vital  to  the  National 
Interest  (MAVNI)  program  this  year.  A  new  student,  Bernard  Moore 
joined  our  group  in  the  spring.  He  participated  the  US  Air  Force  SFFP 
program  and  is  highly  motivated  by  our  project  as  well  as  the  service  to 
the  military.  The  numerical  techniques  we  studied  so  far  have  had  an 
impact  on  education  and  training.  In  the  previous  project,  we  recruited 
and  advised  three  high  school  students  for  their  summer  research  in 
computational  simulation  of  parachute  and  other  fluid-structure  inter¬ 
actions.  These  students  were  Andrew  Ahn  and  Thomas  Hobbs  from 
the  Ward  Melville  High  School,  and  Heisler,  Alicia  from  North  Shore 
Hebrew  Academy  High  School. 
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